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A nonlinear electromagnetic scattering problem is studied in the presence of bound states in the 
radiation continuum. It is shown that the solution is not analytic in the nonlinear susceptibility and 
the conventional perturbation theory fails. A non-perturbative approach is proposed and applied 
to the system of two parallel periodic arrays of dielectric cylinders with a second order nonlinear 
susceptibility. This scattering system is known to have bound states in the radiation continuum. 
In particular, it is demonstrated that, for a wide range of values of the nonlinear susceptibility, the 

j>^ " conversion rate of the incident fundamental harmonic into the second one can be as high as 40% 

when the distance between the arrays is as low as a half of the incident radiation wavelength. The 

f^*) ■ effect is solely attributed to the presence of bound states in the radiation continuum. 
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I. INTRODUCTION 

A conventional approach to nonlinear electromagnetic scattering problems is based on the power series expansion 
in a nonlinear susceptibility Xc- For example, for the 2nd order susceptibility, the physical parameter that determines 
nonlinear effects is XcEr ^ 1 where Er is the electric field at the scattering structure. The smallness of XcEr justifies 
the use of perturbation theory and the solution is analytic in Xc- The situation is different if the scattering structure 
has resonances. 

Planar periodic structures (e.g., gratings) are known to exhibit sharp scattering resonances when illuminated by 
electromagnetic waves (for a review see, e.g., [il0|)- Furthermore, it is known (see, e.g., [2:]) that in such structures 
a local electromagnetic field E^. is amplified if the structure has narrow resonances: E^. ^ Eil\fV , where Ei is the 
amplitude of the incident wave and F is the width of the resonance. Consequently, optical nonlinear effects are 
amplified if the system has a sufficiently narrow resonance. For example, an amplification of the second harmonics 
generation by a single periodic array of dielectric cylinders [3j and by other single-array periodic systems [J] have been 
reported, but no significant fiux conversion rate, comparable to that in conventional methods of second harmonic 
generation, has been found. Owing to the smallness of Xc and a fmite I/a/F, the condition Xc^r ^ 1 still holds for 
the studied structures. 

If two identical planar periodic structures are aligned parallel and separated by a distance 2/i, then it can be shown 
that for each resonance associated with the single structure, the combined structure has two close resonances whose 
width depends continuously on h so that the width of one of the resonances vanishes, i.e., F(/i) -^ as h ^>- hh for 
some discrete set of distances hj,, for sufficiently large h [3,|5|, |6|. This means that the system has bound states in the 
radiation continuum. Their existence was first predicted in quantum mechanics by von Neumann and Wigner [7| in 
1929 and later they were discovered in some atomic systems Q (see also @-[ll| for more theoretical studies). Their 
analog in Maxwell's theory has only attracted attention recently |5l. [l24l4j . In particular, for a system of two parallel 
arrays of periodically positioned subwavelength dielectric cylinders (depicted in Panel (a) of Fig. [1]), the existence of 
bound states in the radiation continuum has been first established in numerical studies of the system ]^ . A complete 
classification of bound states as well as their analytic form for this system is given in [6] for TM polarization. It is also 
shown [y| that bound states exist in the spectral range in which more than one diffraction channel are open. From 
the physical point of view, bound states in the radiation continuum are localized solutions of Maxwell's equations like 
waveguide modes, but in contrast to the latter their spectrum lies in the spectrum of scattering radiation (diffraction) 
modes. 

The perturbation theory parameter XcEr ~ XcEi/\/T{h), as defined above, can no longer be considered small if 
bound states in the radiation continuum are present. This qualitative assessment should be taken with a precaution. 
In the present study, a rigorous analysis of the nonlinear scattering problem by means of the formalism of Siegert 
states (appropriately extended to periodic structures il5|]) shows that no divergence of a local field occurs as /i — )■ /i;,. 
However, the conventional perturbative approach fails because the solution is not analytic in Xc- The situation can 
be compared with a simple mechanical analog. Consider a scattering problem for a particle on a line in a hard core 
repulsive potential V{x) — g/x'^, 3 > 0. No matter how small g is, the particle never crosses the origin x — and a 
full reflection occurs, but it does so when 5 = (a full transmission). So, the scattering amplitude is not analytic in 
g. Other, more sophisticated, examples of quantum systems with such properties are studied in 16]. 

The purpose of the present study is twofold. First, the nonlinear scattering problem is studied in the presence of 
bound states in the radiation continuum. A non-perturbative approach is developed to solve the problem. Second, 
as an application of the developed formalism, the problem of the second harmonic generation is analyzed with an 
example of the system depicted in Fig. [1] (Panel (a)). 



A. An overview of the results 

In Section [TTl the nonlinear resonant scattering problem with bound states in the radiation continuum is trans- 
formed into a system of integral equations. A non-perturbative method is proposed to solve these equations in the 
approximation that takes into account two nonlinear effects: a second harmonic generation in the leading order of Xc 
and the fundamental harmonic generation by mixing the second and fundamental harmonics in the leading order in 
Xc- This is the second order effect in Xc known in the theory as the optical rectification. The latter is shown to be 
necessary to ensure the energy flux conservation. 

The formalism is illustrated with an example of two parallel periodic arrays of dielectric cylinders shown in Panel 
(a) of Fig.[TJ The analysis is based on the subwavelength approximation (Section IIIip when the incident wave length 
is larger than the radius R of the cylinders. If k is the magnitude of the wave vector, then the theory has three small 



parameters: 

So{k) = I {kR)\ec - I) ^ I , Xc«l, \Ah\ = \h - hb\ ^ 1 

where all the distances are measured in units of the structure period, in particular R < 1/2. (5o(fc) is the scattering 
phase for a single cylinder, Sc is the linear dielectric susceptibility, and the amplitude of the incident wave is set to 
one, Ei = 1. With this choice of units, all three parameters are dimensionless. The scattering amplitudes of the 
fundamental and second harmonics are explicitly found in Section [IVI 

In Section [Vj it is shown that the ratio of the flux of the second harmonic along the normal direction and the 
incident flux is 

a2 = CxlEt 

where E^ = Er{xc, A/i, So) is the fundamental field on the cylinders, and C — C{So, Ah) is some function. The function 
Ej, is non-analytic in the vicinity of zero values of its arguments. The non-perturbative approach of Section|lT]is used 
to prove that the generated flux of the second harmonic attains its maximal value when the small parameters satisfy 
the condition 

{Ah)^S',{h)^dxl (1) 

where i? ^ 1 is a numerical constant, and kb is the magnitude of the wave vector of the bound state that occurs at 
h = hfj. Under this condition, (72 becomes analytic in the scattering phase Sq so that in the leading order, 

0-2, max ~ 47rfcf/^(5o(fcfc) 

An interesting feature to note is the independence of the conversion efficiency on the nonlinear susceptibility Xc (in 
the leading order in the scattering phase 5q). In other words, given a nonlinear susceptibility Xc by a fine tuning of 
the distance between the arrays one can always reach the maximal value which is only determined by the scattering 
phase at the wave length of a bound state. The lowest value of kb for the system considered occurs just below the 
first diffraction threshold (the wavelength is slightly larger than the structure period) Q, i.e., kb « 27r. Taking, for 
example, R — 0.15 and £c = 2 (so that 5q{2t:) « 0.22), the conversion rate reads fT2,max ~ 0.44, that is, about 44% of 
the incident flux is converted into the second harmonic flux, which is comparable with the conversion rate achieved 
in slabs (crystals) of optical nonlinear materials [I3| ■ 

From the physical point of view, the scattering structure plays the role of a resonator with the quality factor inversely 
proportional to T. The field in the resonator is not uniform and has periodic peaks of the amplitude Er ~ Eij \fV 
due to a constructive interference of the scattered fundamental harmonic. The second harmonic is produced by the 
induced dipole radiation of point scatters located at these peaks. The induced dipole strength is proportional to 
Xc-E'r- The dipoles are excited by the incident wave and, due to their periodic arrangement, they radiate in phase 
producing a plane wave in the asymptotic region (just like a phased array antenna). If the system has a resonance 
whose width F can be continuously driven to zero by changing a physical parameter of the system, i.e., the system has 
a bound state in the radiation continuum, then the strength of the induced dipoles radiating the second harmonics 
can be magnified as desired, but the resonator cannot be excited by the incident radiation if F = (a bound state 
is decoupled from the radiation continuum). So, the optimal width F at which the second harmonic amplitude is 
maximal occurs for some F ^^ 0, which explains the existence of conditions like ([T]). Since the second harmonic is 
generated by point scatterers, the phase matching condition, needed for optically nonlinear crystals, is not required. 
The energy flux of the incident radiation is automatically redistributed and focused on the scatterers owing to the 
constructive interference. Thanks to these physical features, an active length at which the conversion rate is maximal 
is close to 2/ifc whose smallest value for the system studied is roughly a half of the wave length of the incident light 
p, |6| (i.e. for an infrared incident radiation it is about a few hundreds nanometers). 

II. THE NONLINEAR RESONANT SCATTERING THEORY 

Suppose that a scattering system has a translational symmetry along a particular direction and has non-dispersive 
linear and second-order nonlinear dielectric susceptibilities, e and x^ respectively. When the electric field is parallel 
to the translational symmetry axis (TM Polarization), Maxwell's equations are reduced to the scalar nonlinear wave 
equation 

^9,^(^£^+|^^^)-Ai? (2) 



Let the coordinate system be set so that the functions e — 1 > and X > have support bounded in the z— direction 
and the system has the translational symmetry along the y— direction. In this case, e, x E'-nd E are functions of z 
and X. In the asymptotic regions \z\ — >■ cx), Eq. ([2]) becomes a hnear wave equation. So, the scattering problem 
can be considered for a plane wave of the frequency uj that propagates from the asymptotic region z — ;> — oo to the 
region z — ^ cxd. Furthermore, it is assumed that the functions e and x ai'e piecewise constant, i.e., e ~ Sc ~ const and 
X — Xc — const in regions occupied by the scattering system. A conventional treatment of the problem is based on 
the assumption that the solution E is analytic in Xc and, therefore, can be represented as a power series expansion, 

E = 2Re {Eie-'^' + XcE2C-^"^' + xl (^s.ie"*"* + i^s.ae-^'"*) + . . .} (3) 

where Ei is the amplitude of the fundamental harmonics in the zero order of Xc, E2 is the amplitude of the second 
harmonics in the first order of Xc and so on. This assumption is not true if the system has bound states in the 
radiation continuum. Indeed, a general solution has the form E — E^ + E^l, where E^ is the solution when Xc = 
and Enl is the correction due to nonlinear effects. Let x be written as x = Xc^j where 77 is the indicator function 
of the region occupied by the scattering system, i.e., its value is 1 in that region and elsewhere. Then, if G is 
the Green's function of the operator ■^d'f — A with appropriate (scattering) boundary conditions, the function Eiq^ 
satisfies the integral equation 

Enl = -^G [vd^tiEL + Enl?] 

The power series expansion ([3|) can be obtained by the method of successive approximations for this integral equation, 
provided the series is proved to converge. According to scattering theory [1^, [l^, the Fourier transform of G is 
meromorphic in A:^ = ^. As is clarified shortly (see discussion of Ea. ([TT|) ). its real poles correspond to bound 

states in the radiation continuum. Hence, in the presence of a real pole k^ — fc^, the kernel of G is not summable 
and, therefore, the successive approximations produce a diverging series. This implies a non-analytic behavior of the 
solution in Xc- Thus, when a bound state in the radiation continuum is present, the conventional perturbative approach 
becomes inapplicable. Here, a non-perturbative approach is developed to obtain the solution to the scattering problem 
that is valid in any small neighborhood of a real pole of the Fourier transform of G. 
Suppose that the incident radiation is a plane wave 

i?i„(r,i) = 2cos(k • r — wi), k = fc^jGi + ^263, ck = Ld, 

where e.^, i — \,2, 3, denote unit vectors along the x, j/, and z coordinate axes, respectively. A general solution to Eq. 
^ should then be of the form, 

00 
E{v,t)= Y, Ei{v)c-'"^' (4) 

/— — 00 

where Eq = 0, and for all I, E-i — Ei is the complex conjugate of Ei (as E is real). Therefore it is sufficient to 
determine only Ei, I > 1. Next, it is assumed that the scattering structure is periodic in the a;— direction (e.g., a 
grating). The units of length are chosen so that the period is one. Then the amplitudes Ei satisfy Bloch's periodicity 
condition 

Sz(r + ei) = e^"==Sz(r) (5) 

This condition follows from the requirement that the solution E satisfies the same periodicity condition as the incident 
wave Ein'- 



ijj 
By Eq.©, the amplitudes of the different harmonics satisfy the equations, 



£'„(r + ei,i) =Ein\v,t ^j 



l\Ei + l^k^eEi = -vl^k^{e - 1) ^1 ^p^i-p^ 



M^c - 1) 



For ease of notation, the parameter ly is often used in lieu of Xc- Since t/ ^ Xci it is a small parameter in the system. 

The scattering theory requires that for I ^ ±1, the partial waves Eie~^^'^* be outgoing in the spatial infinity 

(|z| — >■ 00). The fundamental waves iSj-ie^*'^* are a superposition of an incident plane wave ^^^i^-'^-'^t) and a 



scattered wave which is outgoing at the spatial infinity. In all, the above boundary conditions lead to a system of 
Lippmann-Schwinger integral equations for the amplitudes Ei: 



(6) 



- B,{{lky)[Ei+:yj:^EpEi^p], l>2 
and E-i = Ei ioi I < —1, where H(g^) is the integral operator defined by the relation 

H(g2)[^](r) = ^ J{e{ro) - l)G,(r|ro)V(ro)dro (7) 

in which G^(r|ro) is the Green's function of the Poisson operator, {q^ + A)Gq(r|ro) = — 47r(5(r — Tq), with the outgoing 
wave boundary conditions. For two spatial dimensions, as in the case considered here r = {x, z) and ro = (cco, zq), the 
Green's function is known [20] to be Gg(r|ro) = iTTHf){q\Y — ro|) where Hq is the zero order Hankel function of the 
first kind. 

When i^ = 0, the amplitudes of all higher harmonics (/ > 2) vanish. Therefore it is natural to assume that 
\Ei\ ^ \E2\ ^ I-E3I ^ • • • for a small v. Note that this does not generally imply that the solution, as a function of 
I/, is analytic ai u = Q. Under this assumption, the solution to the system ^ can be approximated by keeping only 
the leading terms in each of the series involved. In particular, the first equation in (|5]) is reduced to 

El « e^'"'- + 'R{k')[Ei\ + 2!.H(fc2) [E1E2] (8) 

while the second equation becomes 

E2 « il((2fc)2)[£;2] + z.H((2fc)2) [El] (9) 

It then follows that a first order approximation to the solution of the nonlinear wave equation (^ may be found by 
solving the system formed by the equations ([8]) and ([9]). To facilitate the subsequent analysis, the system is rewritten 
as 

f[l-H(fc2)][iJ,]=e*k-"- + 2z.H(fc2)[Eii?2] , . 

\[l-m2kn[E2] = uYi{{2kf)[El] ^ ^ 

Solving the first of Egs. pUl) involves inverting the operator 1 — Il(fc^), and therefore necessitates a study of the 
poles of the resolvent [1 — H(fc^)]~^. Such poles are eigenvalues in the generalized eigenvalue problem, 

n{k^)[E]^E (11) 

for fixed k^. The corresponding eigenfunctions E — Eg are referred to as Siegert states. In contrast to Siegert states 
in quantum scattering theory jl8| , electromagnetic Siegert states satisfy the generalized eigenvalue problem (II ip in 
which the operator is a nonlinear function of the spectral parameter k^ . Their general properties are studied in [l^ . 

Eigenvalues have the form k"^ = k^ — iV. If kr > k^, then, according to scattering theory, such a pole is a resonance 
pole. In the case of the linear wave equation {xc = ^ = 0), the scattered flux peaks at k = kr indicating the resonance 
position, whereas the imaginary part of the pole F defines the corresponding resonance width (or a spectral width of 
the scattered flux peak; a smaU F corresponds to a narrow peak). If F = 0, the corresponding Siegert state is a bound 
state. This is a localized (square integrable) solution of Eq. pTj) . Suppose that the scattering system has a physical 
parameter h such that a pole k'^ = k'^{h) — iT{h) of the resolvent [1 — H(k'^)]~^ depends continuously on h and that 
there is a particular value h = hb at which the pole becomes real, i.e., T{hi,) = 0. If kb = kr{hb) > k^, then the 
corresponding bound state lies in the radiation continuum. Note that Eq. lTTI) may have solutions for real k^ < k^. 
These are bound states below the radiation continuum. Such states are not relevant for the present study and, 
henceforth, bound states are understood as bound states in the radiation continuum. As noted in the introduction, 
two periodic planar scattering structures separated by a distance 2h have bound states in the radiation continuum. 

Suppose that for fixed k^, the set k^ consists of isolated points (which is generally true) and the points (2fcf,)^ do 
not belong to it (which is fulfilled in a concrete example studied in the next section G]). Consequently, the operator 
(1 — H((2fc)^)~^ is regular in a neighborhood of k^ and, for k close to but not equal to kb, the operator 1 — H(A:^) is 
invertible. It then follows that Ei satisfies the nonlinear integral equation. 



El = (l-H(fc2 



e"^"- + 2t/^H(fc^) 



Ei(i~Rii2kr)y' R{{2kr)[E-t\ 



(12) 



where, in accord with the notation introduced in dT]), the function on which an operator acts is placed in the square 
brackets following the operator. The operator iv^(l — H(A:^))~^ that determines the "nonlinear" part of Eq. ([T2|) is 
not bounded as k^ —> fc^, no matter how small i/^ ^ xi is- This precludes the use of a power series representation 
of the solution in i/^. To find the solution of Eq. (ll2|) when ly is small, its property under parity transformations is 
established first. ^ ^ ^ 

Suppose that the scattering system is such that the operator H and the parity operator P defined by P[i?](x, z) — 
E(x, —z) commute. This implies that the Siegert states have a specific parity: P[i?s] = pEs where p = ±1. Consider 
then the ratio 

,i...,.Ei^ fm (13, 

Ei{x,+z) El 

It will be proved that /i(x, z) — > p in the limit {h,k) — > {hb,ki,) along a certain curve. Indeed, it follows from the 
meromorphic expansion of [1 — H(fc^)]^^ that near a pole k^{h) — ir{h), 

where C{h) is some constant depending on h, and Eg is an appropriately normalized Siegert state [l5|. Consider then 
the curve of resonances C : k — kr{h) in the {h, fc)-plane. Along C, 

Ei{x,z)^^E,{x,z) + 0{l) (15) 

Now, as ft, — > hi,, the width T{h) goes to 0, and the Siegert state Eg becomes a bound state Ei, in the radiation 
continuum. Equation p^ shows that if C{h) does not go to zero faster than r(/i) as ft, — >■ ftf,, i-e., the pole still gives 
the leading contribution to Ei in this limit, then 

Eb{x,-z) P[Eb] 
Kx,z)^ =—-^=p = ±l 16 

Eb{x,+z) Eb 

depending on whether the bound state Eb is even or odd in z. For the linear wave equation (j/ == 0), the constant 
C(ft) is shown to be proportional to ^y^(h) [15|. Therefore, for a small v, the assumption that C{h) does not go to 
zero faster than r(ft) as ft — > ftb is justified. 

Based on the limit (|16p . the following (non-perturbative) approach is adopted to solve Eg. p^ near a bound state. 
First, the curve of resonances C in the (ft, fc)-plane is found. Using the relation P^ in the right side of P^ . the field 
El is expressed via the ratio /z with the pair (ft, k) being on the curve C. Next, the principal part of the amplitude 
El relative to Aft = ft — ftf, is evaluated near a critical point {hb,kb) on 6 by taking /i to its limit value ([TO)) . This 
approach reveals a non-analytic dependence of the amplitude Ei on the small parameters Aft and Xc of the system 
and allows to obtain Ei and E2 when a bound state is present in the radiation continuum. As the technicalities of 
the proposed non-perturbative approach depend heavily on peculiarities of the scattering system, the procedure is 
illustrated with a specific example. 

III. A PERIODIC DOUBLE ARRAY OF SUBWAVELENGTH CYLINDERS 

The system considered is sketched in Fig. [TJa). It consists of an infinite double array of parallel, periodically 
positioned cylinders. The cylinders are made of a nonlinear dielectric material with a linear dielectric constant £c > 1, 
and a second order susceptibility Xc ^ 1- The coordinate system is set so that the cylinders are parallel to the y-axis, 
the structure is periodic along the x-axis, and the z-axis is normal to the structure. The unit of length is taken to be 
the array period, and the distance between the two arrays relative to the period is 2ft. 

The solution of the integral equation ((T2|) is obtained for fe^ near k^ in the limit of subwavelength dielectric cylinders. 
The approximation is defined by a small parameter 

^o(g) = ^(£c-i)«l (17) 

which is the scattering phase of a plane wave with the wavenumber q on a single cylinder of radius R. For sufficiently 
small R, this approximation is justified. The integral kernel of Il(q^) is defined by (O and has support on the region 
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FIG. 1. Panel (a): Double array of dielectric cylinders. The unit of length is the array period. The axis of each cylinder is 
parallel to the y-axis, and is at a distance h from the x-axis. 

Panel (b): The scattering process for the normal incident radiation {kx = 0). The scattered fundamental harmonic is symbolized 
by a single headed arrow while the (generated) second harmonic radiation is symbolized by a double headed arrow. The incident 
radiation wave length is such that only one diffraction channel is open for the fundamental harmonic while three diffraction 
channels are open for the second harmonic. The flux measured through the faces Z/±i/2 : x = ±1/2 cancels out due to the 
Bloch periodicity condition as explained in Appendix [B] 

Panel (c): The solid and dashed curves show the position (frequency Ur = cfc,) of scattering resonances as functions of the 
distance between the arrays, k = kr{h). The dots on the curves indicate positions of bound states in the radiation continuum 
(i.e., the values of h at which a resonance turns into a bound state). The solid curve connects bound states symmetric relative to 
the reflection z — )■ —z. The dashed line connects the skew symmetric bound states. The curves are realized for R — 0.08, Sc ~ 2, 
and kx = (normal incidence). 

occupied by cylinders. The condition ()17p implies that the wavelength is much larger than the radius R, and therefore 
field variations within each cylinder may be neglected, so that ip{x, z) ~ i^{n, ±h) where (n, ±/i) are the positions of 
the axes of the cylinders {n is an integer). The integration in Ji{q'^)[ip] yields then an infinite sum over positions of 
the cylinders. By Bloch's condition, i/j{n,zLh) = e™''^7/;(0,±/i), so that the function ii{q'^)[ip]{x, z) is fully determined 
by the two values ip{0, ±h). In particular. 



H(g2)[^](0, ±h) « aV'(0, ±h) + ^V(0, Th) 



(18) 



where the coefficients a and /3 are shown to be 



a{q,qx) = 2TTiSa{q) 




2Tri{\m\ 



- ln(27ri?) 



I3{q,qx,h) ^2TTi6o{q) ^ 



771 — — OO 



Qz,m 



where q^^m — \/q^ — {Qx + ^irmY with the convention that if q^ < [qx + 27rm)^, then qz^m — i\/[qx + 27rm)2 — q^. 
To obtain the energy flux scattered by the structure, the action of the operator H((7^) on -0 must be determined in 
the asymptotic region |z| -^ oo. It is found that for \z\ > h + R, 



OO 

U{q^M{x, z) « 2^i6o{q) ^ (^'(O, /i)e*l--''l'?-'" + 7A(0, -/i)e'l^+''l«-'") 



771 — — OO 



oia:((j^+27rm) 



(20) 



IV. AMPLITUDES OF THE FUNDAMENTAL AND SECOND HARMONICS 



Now that the action of the operator H(g^) has been established in p^ and (^0)) . the amplitudes Ei and E2 of the 



fundamental and second harmonics can be determined by solving the system pH]) . As noted earlier, this will be done 
along a curve C in the /i, A:-plane defined by A: = kr{h) where kr{h) is the real part of a pole of [1 — H(fc^)]~^, or 
equivalently, when the incident radiation has the resonant wave number k = kr{h). To find the curve, the eigenvalue 



problem pT|) is solved in in the approximation ([18 



where Eb± ~ Eb{0,±h) and the functions a — a{k,kx) and /3 = f5{k,kx) have been defined in the previous section. 
In particular, bound states occur at the points {hi,, k^) at which the determinant det(l — 3^) vanishes, 



det {^_^ ^ M = (1 - a - /?)(! - a + /3) = 



It follows from Eq.(|2T|) that the bound states for which 1 — a — /? = are even in z because Eb+ = Eh- in this 
case. Similarly, the bound states for which 1 — a + /3 = are odd in z. More generally, the poles of the resolvent 
[1 — H(A;^)]~'^ are complex zeros of det(l — 'K). They are found by the conventional scattering theory formalism. 
Specifically, the resonant wave numbers k^ — k^{h) are obtained by solving the equation Re{l — a ± /3} = for the 
spectral parameter fc^. According to the convention adopted in the representation (1141) . the corresponding resonance 
widths are defined by 



r(,)^_.Mi-"±/3} 



dk2Re{l-a±P} 



k'2=kl{h) 



where 9^2 denotes the derivative with respect to fc^. This definition of the width corresponds to the linearization 
of Re{l — a ± /3} near fc^ = fc^(/i) as a function of k^ in the pole factor [1 — a ± /3]~^. The curves of resonances 
k = kr{h) > kx come in pairs. There is a curve connecting the symmetric bound states in the h, fc-plane, and another 
curve that connects the odd ones. 

In what follows, only the curve connecting symmetric bound states will be considered. The other curve can be 
treated similarly. Panel (c) of Fig. [T] shows that the first symmetric bound state occurs when the distance 2h is 
about half the array period, while the skew-symmetric bound states emerge only at larger distances. This feature is 
explained in detail in 6]. So, the solution obtained near the first symmetric bound state corresponds to the smallest 
possible transverse dimension of the system (roughly a half of the wave length of the incident radiation). Thus, from 
now on the curve of resonances 6 refers to the curve in the /i, fc-plane defined by the equation Re{l — a — /3} = 0. To 
simplify the technicalities, it will be further assumed that only one diffraction channel is open for the fundamental 
harmonics, i.e., kx < k <2tt — kx- 

Let k — kr{h) be the solution of Re{l — a — /3} = 0. By making use of the explicit form of the functions a and /3 
for one open diffraction channel, one infers that along the curve k = kr{h), 



1 — a — /? = ilm{l — a ~ /3} — — i ip , ip — cos{hkz), kz = y k^ — fc^ 



Bound states in the radiation continuum occur when the distance h satisfies the equation ip ~ 0, i.e., h^Jk'^{K) — k^ — 
{n — l/2)7r with n being a positive integer. Its solutions h ~ hi,{n) define the corresponding values of the wave numbers 
of the bound states, kh{n) — kr{hi,{n)). So, the sequence of pairs {(/ih(n), fc(,(n))}^]^ indicates positions of the bound 
states on the curve 6. In the limit h — > hi,{n) along C, the function ip has the asymptotic behavior, 

(^ =: (-l)"fc^,f,A/i + o(A/i), kz^b = kz\h=ht(n}, Ah^h-hb 

The objective is to determine the dependence of the amplitudes Ei and £'2 on the parameters Ah and Xc which are 
both small. 

To this end, let Ei± = Ei{0,zLh) be the values of the field Ei on the axes of the cylinders at (0,±/i), and 
E2± — i?2(0,±/i) be the values of the field E2 on the same cylinders. In the subwavelength approximation, these 
values determine the scattered field because the latter is produced by the radiation of point dipoles induced by the 
incident wave on the scatterers and the strength of the dipoles is proportional to Ei± for the fundamental harmonic 
and i?2± for the second harmonic. Applying the rule ([T51) to evaluate the action of the operator Il(g^) in the system 
PH)) . the first equation of the latter becomes, 

ii-Ki(M^2.«(|»M+(^r,;::) (22., 



Similarly, the second of Eqs. dTO)) yields. 



[1 - K2] (fA = ^^^2 (^+ ], ^2= n^k, 2kx) (22b) 



As stated above, the resolvent [1 — H((2fc)^)]~^ is regular in a neighborhood of ki, so that Eq. (j22bl) can be solved 
for i?2±, which defines the latter as functions of Ei±. The substitution of this solution into Ea. (|22ap gives a system 
of two nonlinear equations for the fields Ei±. Adding these equations and replacing the field Ei- by its expression 
El- = /i(0, h)Ei+ in terms of the field ratio of Eq.((T3]) yields the following implicit relation between the field Ei+ 
and its amplitude |i5i+|: 

Ei+ = 2 (23) 

where Lp — cos{hkz) and v — ^^,'^''_j^x are small and, in terms of the field ratio /i = //(O, h), the values of C and ^ read, 



with a and b being defined by the relation. 



b 



1-:K,]-':K,^{1 ;) (25) 



In particular, ^ and ^ are continuous functions of ^ and (p. In Appendix|X]it is shown that if ^6 and £,b are the respective 
limits of C and ^ as /i — )■ /i;, along the curve of resonances C, then these limits are nonzero. It follows then that Eg. ([25)) 
for El-^- is singular in both v and ip when these parameters are small, i.e., in the limit (j^, <p) — > (0,0). Furthermore, 
there is no way to solve the said equation perturbatively in either of the parameters. A full non-perturbative solution 
can be obtained using Cardano's method for solving cubic polynomials. Indeed, by taking the modulus squared of 
both sides of the equation, it is found that, 

2 4 2 

X^ + 2^Re{a}X^ + ^\CeX~^\C\'^0, X^\Ei+\' (26) 

The solution to this cubic equation is obtained in Appendix [C] It is proved there that Eq. ([26]) admits a unique 
real solution so that there is no ambiguity on the choice of -Bi+. In the vicinity of a point {hi,, kh) along the resonance 
curve 6, the field i?i+ is found to behave as, 

\Ei+\ = L^r{Ah,Xc) (27) 

Xc 

Recall that Ah = h — hi,. An explicit form of the function r(Aft,, Xc) is given in Appendix [Cj f see Eq. (|C3p ). It involves 
combinations of the square and cube roots of functions in Ah and Xc and has the property that t(A/i, Xc) — > as 
{Ah,Xc) -^ (0,0) (in the sense of the two-dimensional limit). In the limit h — > hf,, the field ratio /i approaches 1 for 
a symmetric bound state as argued earlier. Therefore it follows from Eqs. (l22b[) that E2± ~ i^^i+ because the matrix 
(I^SI) exists at h — hf,. Since v ^ Xc, relation (P7|) leads to the conclusion that 

fi.O(|A/.|V3) 

Ei+ 



Thus, the approximation \Ei\ ^ \E2\ used to truncate the system (fTO| remains valid for h close to the critical value 
hf, despite the non-analyticity of the amplitudes at {Ah,Xc) ~ (0,0). 

V. FLUX ANALYSIS: THE CONVERSION EFFICIENCY 

For the nonlinear system considered, even though Poynting's theorem takes a slightly different form as compared 
to linear Maxwell's equations, the flux conservation for the time averaged Poynting vector holds. The scattered 
energy flux carried across a closed surface by each of the different harmonics adds up to the incident flux across that 
surface. The flux conservation theorem is stated in Appendix [B] Consider a closed surface that consists of four faces, 
L± = {{x,z)\ — ^ < X < \,z ^ ±00} and i±i/2 = {{x,z)\x — ±1/2} as depicted in Fig.[Ub). As argued in Appendix 
IB| the scattered flux of each Z'^'-harmonics across the union of the faces L±i/2 vanishes because of the Bloch condition 
(and so does the incident flux for any k^). Therefore only the flux conservation across the union of the faces L± has 
to be analyzed. If ct; designates the ratio of the scattered flux carried by the /"^-harmonics across the faces L± to the 
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incident flux across the same faces, then J2i>i '^i — 1- Thus, for I > 2, ai defines the conversion ratio of fundamental 
harmonics into the Z*''-harmonics. 

In the perturbation theory used here, only the ratios ci and (72 niay be evaluated. By laborious calculations it can 
be shown that ci + (T2 < 1 as one would expect (see Appendix IB] for details). Hence, the efficiency of converting the 
fundamental harmonic into the second harmonic is simply determined by the maximum value of a2 as a function of 
the parameter /i at a given value of the nonlinear susceptibility Xc- 

The ratio (72 is defined in terms of the scattering amplitudes of the second harmonic, i.e., by the amplitude of E2 
in the asymptotic region \z\ — > 00: 



^2(r) 






'*''e*'■■•'".=^ z- 



s/i ir-k^ 



z — > +00 



(28) 



where k^ 



{2kx + 27rm)ei ± fc|^e3 is the wave vector of the second harmonic in the m open diffraction channel. 



Recall that the m**^ channel is open provided (2fc)^ > {2kx + 27rm)^ and in this case fc|''^ = \J{2kY 
while if the channel is closed, then fcf'j^ = iyj{2kx + 27r?7i)2 — [2kY. In the asymptotic region \z\ — 



(2fca; + 27rm)2, 
cxD, the field in 



iyj{2kx + 27rTO)2 

closed channels decays exponentially and, hence, the energy fiux can only be carried in open channels to the spatial 
infinity. The summation in Eas. (j28p is taken only over those values of m for which the corresponding diffraction 
channel is open for the second harmonic, which is indicated by the superscript "op, s/i" in the summation index 
m"^'"^ . Note that there is more than one open diffraction channel for the second harmonic even though only one 
diffraction channel is open for the fundamental one. For instance, if the x— component of the wave vector k, i.e., 
kx^ is less than ^, there are 3 open diffraction channels for the second harmonic, the channels m = 0, m = —1, and 
m = 1. These three directions of the wave vector of the second harmonic propagating in each of the asymptotic regions 
z — ?> ±cx) are depicted in Fig. [IJb) by double-arrow rays. Thus, in terms of the scattering amplitudes introduced in 
Eqs. ([28)) . the ratio of the second harmonic flux across L± to the incident flux is 



''^-^ T. K%A\R:n\' + \T, 



sh\2\ 



The scattering amplitudes RfJ^ and T^ are inferred from Eq.® in which the rule (PO]) is applied to calculate the 
action of the operator II((2A:)^) in the far-field regions \z\ —>■ cx): 



R 



,,, _ 2Tn6oi2k) 



Ush 



j,sh ^ 2mSo{2k) 



Ush 



{E2+ + vEl^)e' 



{E2+ + ^.i?i\)e-""=^'™ 



{E2- + vEl_)e-''''''^- 
+ (S2-+i^^L)e'''^- 



Since for v ^ Q, the amplitudes Ei± remain finite as /i — > /ib along C, and since /i — )■ 1 in the said limit; the principal 
part of 0-2 in a vicinity of a bound state along the curve C is obtained by setting £'i_ — _Ei+ in Eq. (j22bp . solving the 
latter for E2±, and substituting the solution into the expression for cr2- The result reads 



a2 ^ Ci,v^\E^+\'' 

where Cb is a constant obtained by taking all nonsingular factors in the expression of a2 to their limit as h 
which gives 



(29) 



Cb 



(16^(5o(fc))^ 



b\ 



^E 

^op.sh 



cos^ (hk 



sh \ 
z.m) 



Ush 
'^z.m 



{h,k) = {hi,,kb) 



for a and b defined in Eq. (P5|) . Using the identity |i?i+|^ = |i?i+p|i?i+p, and substituting Eq. (P5)) into one of the 
factors l-Bi+p, the conversion ratio a2 is expressed as a function of a single real variable, 

<J2iu) = C'b- 



l^ + Cfc^P' 



iy\Ei^ 



f 



(30) 



CblCbP is a constant, and ^ and C in Eq. ([M)) have been taken at their limits a.s h —> hb to obtain the 



where C^ 

principal part of (T2- The function u 1— > o'2{u) on [0,oo) is found to attain its absolute maximum at u — |^bCb|- This 
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condition determines the distance 2h between the arrays at which the conversion rate is maximal for given parameters 
R, Ec and Xc of the system. Indeed, since u ~ |^i+P should also satisfy the cubic equation (^5]) . the substitution of 
u = IC&Cbl iiito the latter yields the condition 



— = 2|6P(l6ai+Re{6Cfc}) 



(31a) 



In particular, in the leading order in <5o(^), the optimal distance 2h between the two arrays is given by the formula, 

,2 



(h - h„Y 



Xc 



?7r5A:,,b(fcbi?)6(e,-l)5 



(31b) 



where as previously, k^^ = \f% 






The maximum value (72, max of the conversion ratio (72 is the sought-for conversion efficiency. An interesting feature 
to note is that a2,max = cr2{\Cb^b\) is independent of the nonlinear susceptibility Xc because the constants C^, ^6, 
and Cb are fully determined by the position of the bound state (fcb, /if,). In other words, if the distance between the 
arrays is chosen to satisfy the condition (I31ap . the conversion efficiency (72, max is the same for a wide range of values 
of the nonlinear susceptibility Xc- This conclusion follows from two assumptions made in the analysis. First, the 
subwavelength approximation should be valid for both the fundamental and second harmonics, i.e., the radius of 
cylinders should be small enough. Second, the values oi h — hf, and v (or Xc) rnust be such that the analysis of the 
existence and uniqueness of l-Ei+l given in Appendix [Q holds, that is, Eq. (j26p should have a unique real solution 
under the condition pial) . The geometrical and physical parameters of the studied system can always be chosen to 
justify these two assumptions as illustrated in Fig. [5] 




i- -ts 




FIG. 2. Panel (a): The conversion efficiency is plotted against the cyhnder radius R for the critical points {hb{n), kb(n)), n = 

1, 2, 3, and Ec = 1-5 and fc^ = (the normal incidence). The dashed parts of the curves indicate the regions where So{k) > 0.25 

and, hence, 5o{2k) > 1, i.e., the subwavelength approximation becomes inapplicable for the second harmonic. 

Panel (b): The conversion efficiency is plotted against k^ for the critical points {hb{n),kb{n)), n = 1,2,3. The curves are 

realized for R = 0.15 and Sc = 1.5. 

Panel (c): The region of validity of the developed theory for the first bound state hb{l) ~ 0.259. The shadowed part of the 

{v, /i)— plane is defined by the condition t+ + t- > under which, according to Ea. (|C3|) . the amplitude |-Ei+| exists and unique 

as explained in Appendix[Cl The plot is realized for R — 0.1, £c = 2 and k^ = 0. The parabola-like curve is an actual boundary 

of the shadowed region; the top horizontal line represents no restriction. So there is a wide range of the physical and geometrical 

parameters within the shadowed region which satisfy H31ap . The regions of validity for other bound states looks similar. 



Panels (a) and (b) of Fig. [2] show the conversion efficiency f72,max for the first three symmetric bound states 
n = 1, 2, 3, as, respectively, a function of the cylinder radius R when kx = and of kx when R = 0.15. For all curves 
presented in the panels, Sc = 1.5. The values of (72, max are evaluated numerically by Eq. pop where u = \£,bCb\- The 
solid parts of the curves in Panel (a) correspond to the scattering phase do{k) < 0.25 with k — k^. Note that the 
wavelength at which the second harmonic generation is most efficient is the resonant wavelength defined by fc = kr{h) 
where h satisfies the condition (I31ap . For a small Xci the scattering phase at the resonant wavelength can well 
be approximated as i5o(fc) sa (5o(fcf,). The condition (5o(A:fc) < 0.25 ensures that the scattering phase for the second 
harmonic satisfies the inequality (5o(2fcf,) = A5a{ki,) < 1 otherwise the validity of the subwavelength approximation 
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cannot be justified. The dashed parts of the curves in Panels (a) of Fig.[2]correspond to the region where So{ki,) > 0.25. 
Panels (a) and (b) of the figure show that the conversion efhciency can be as high as 40% for a wide range of the 
incident angles and values of the cylinder radius. Such a conversion efficiency is comparable with that achieved in 
optically nonlinear crystals at a typical beam propagation length (active length) of a few centimeters, whereas here 
the transverse dimension 2h of the system studied here can be as low as a half of the wavelength, i.e., for an infrared 
incident radiation, 2h is about a few hundred nanometers. Indeed, as one can see in Fig. [TJc), the first bound state 
occurs at hi,{l) w 0.259 and kf, sa 27r which corresponds to the wavelength Xb = 27r/fcf, w 1. 
The stated conversion efficiency can be fairly well estimated in the leading order of 6o{k): 



0'2,max = a'2(|Cb6l) ~ 8TrSa{k) 



^ cos^hkt^) 



/ J h,sh 



(32) 
(h,k) = {h,,,kf,) 



Suppose that only one diffraction channel is open for the incident radiation. Then m — 0,±1 in Eq. ((32|) (three open 
channels for the second harmonics). Let (tS^„^ denote the term ra = Q, i.e. a2„,„^ is the second harmonic flux in 
which the contribution of the channels with m — ±1 is omitted. In particular, cF2max "^ <''2.max- One infers from Eq. 
(15a that, 

'^2,max ~ , CQs\2hhk,^b) 

As the pair (hb, kb) at which a symmetric bound state is formed satisfies the equation cos{hbkz^b) ~ 0, it follows that 
cos{2hbkz,b) — —1- Hence, 

fc2 

The wavenumbers kb at which the bound states occur lie just below the diffraction threshold 2TT — kx, i.e., kb ~ 2ii — kx 
(see Fig. [TJc) and p\ for details). So that in the case of normal incidence [k^. ~ 0), the above estimate becomes, 

4max«2^M')(£c-l) 

with 5q{2tt) <C 1. If, for instance, R = 0.15 and £c — 2, then, 

^2,nax « 44% 

for 5o{2tt) k. 0.22. 

Appendices 

Appendix A: Estimation of C^ and ^ 

Here the limit values C,b and ^h of the functions C, and ^ defined in Eq. ([M| are estimated as h ^ hb along the 
resonance curve 6. For ^h this is immediate. Indeed, in the aforementioned limit, the field ratio /i — > 1 for a symmetric 
bound state and, therefore, 

. 47rJo(fc) 
th = I — -. 



For C,b, the estimate follows from that wave numbers kb at which bound states exist are close to the diffraction threshold 
2tt — kx when only one diffraction channel is open for the fundamental harmonic, i.e., k^ < k < 27t — k^ [6|. Indeed, 
in the first order of 6o{k), 

ZTT — Kx 

This proximity of the wavenumbers kb to the diffraction threshold 2tt — kx allows one to determine the leading 
terms in the coefficients a and b defined in Ea. (|25l) . and hence Cfo- To proceed, the coefficients a2 = a(2k, 2kx) and 
/32 = I3{2k, 2kx, h) of the matrix IK2 are rewritten by separating explicitly the real and imaginary parts: 

a2 + 132 ^ Tp+ + iSc, a2 - P2 ^ ip- + iS, (A2a) 



13 



where Sc and Sg are defined by the relations, 



S. = 16.Jo(fc) E L ' ^^ - 16.^o(fc) E L ' (A2b) 






and the index rn°P'^^ indicates that the summations are to be taken over all open diffraction channels for the second 
harmonic. Using the estimate (jAip . the functions '0± are found to obey the estimates, 

V'+ = 2 + 0((5o(fcb)), ^_ - 0((5o(fcb)) 

These expressions are then used to estimate C^^^ — 2(a + h). In the first order of So{kb) one infers that 

cos^hk^t,) 



Cb^-\-2m6o{k) J2 



h.sh 



(A3) 

(h,k) = {hb,ki,) 



Appendix B: Complements on the flux analysis: Flux conservation 

For the nonlinear wave equation ([2]), the Poynting Theorem becomes, 



87rdt 



\ 37r 



dv 



S • dn (Bl) 



where T^ is a closed region, and dV is its boundary. The vector S = E x B is the Poynting vector (for simplicity, it is 
assumed that dV lies in the vacuum so that fi = e ~ 1, and x = in a small neighborhood of dV). In the case of a 
monochromatic incident wave, the flux measured is the time-average of S over a time interval T — > oo. By averaging 
Eq.dHH), it then follows that. 



1 '-^ 



(S) ■dn = 0, (S) = - / S{t)dt, T-^oo 

dV ^ Jo 

This is the flux conservation. In terms of the different harmonics of Eq. Q , the time averaged Poynting vector becomes. 

Of interest is the flux of the Poynting vector across the rectangle depicted in Fig. [ijb). By Bloch's condition ([5]), the 
contributions to the flux from the faces L±i/2 : x = ±\ cancel out so that the flux measured is through the vertical 
faces L± = {{x, z)\ — \ < x < ^, z — >■ ±oo}. Note that the vanishing of the flux across the faces i±i/2 is a consequence 
of the fact that the incident wave is uniformly extended over the whole a;— axis. For example, consider the normal 
incidence {k^ — 0) with one diffraction channel open for the incident radiation. Then the Poynting vector of the 
reflected and transmitted fundamental harmonic is normal to the structure and, hence, carries no flux across -/j±i/2- 
The second harmonic {I = 2) has three forward and three backward scattering channels open, m — 0, ±1, relative to 
the z— axis. The wave with m = Q propagates in the direction normal to the structure and does not contribute to 
the flux across L±i/2- Since the incident wave has an inflnite front along the a:— axis, so do the scattered waves with 
TO = ±1. The waves with m = 1 and to, = — 1 carry opposite fluxes across each of the faces i±i/2 as the corresponding 
wave vectors have the same z— components and opposite a;— components and, hence, the total flux vanishes. For a 
finite wave front (but much larger than the structure period), the second harmonic would carry the energy fiux in all 
the directions parallel to the corresponding wave vectors in each open diffraction channel. 

If Gi is as defined in Section|Vl then the flux conservation implies that X^z^i '^i = 1- Therefore, in the perturbation 
theory used, i.e., when the system ([5]) is truncated to Eas. (|10p . the inequality cti + (T2 < 1 must be verified to justify 
the validity of the theory. 

The conversion ratio (72 is given in Section [V] If only one diffraction channel is open for the fundamental harmonic, 
then the ratio cti of the scattered and incident fiuxes of the fundamental harmonic reads, 

ai = |l + roP + |i?o|' 
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where Tq and Rq are the transmission and reflection coefScients which are obtained from the far-field amplitude of 

El as, 



El 



(l+^o)e-^ 



Z — >■ — CXD 

z — > +00 



where k = k^ei + kze^, is the incident wave vector and k — k^ei — kze^ is the wave vector of the reflected fundamental 
harmonic. It then follows from Egs.® and ([20l) that 



Ro^i 

To = i 



2TrSo{k) 



2ii5o{k) 



(Ei+ + 2vE2+Ei+)e'^''' + {E^^ + 2uE2-Ei^)e-'^''' 
{Ei+ + 2uE2+Ei+)e-'^^' + (Si_ + 2vE2-Ei^)e''''^' 

In the vicinity of a critical point (/it,, fc^), the coefficients i?o and To obey the estimate. 



i?o ~ To w i 



4^5o(fcfc) _r. (. . ^''\Ei " 



'<'Z,b 



-V>Ei+ 1 



After some algebraic manipulations, it is found that, 
o"! + 0-2 = 1 



SnSoikh) 2,r^ ,4 
-f \Ei+\ 



fc. 



Ab + , 'V Re _ 



^1+1 



lai 



where At, is the constant defined as, 



Ah 



2Sc\a + b+l\^ -Im 



C ) \ ihM)=(htM) 



and Sc — Im{a2 + P2} is introduced in Eas. (jA2[) . Expressing a and b defined by ((25l) via the coefficients a2 and /32 
of the symmetric matrix 3^2 , one also obtains 



Sc = Im 



a + b 



1 + a + b^ 
Since at the point {hi,, kt) the value of C is Cb = (2(a + 6))^^|(h,fe)=(;it,fe^), it follows that, 

a + b 



A = 



2Im 



l + a + 5 



5+ir-2Im{a + 6} 



{h,k) = {hh,k,,) 

For general complex numbers a and 6, the expression in square brackets is always zero. Therefore Ah = 0, and, 

ai+a2^1 + 2{ , ^ "V t^l-gi+l I I Re' ' ' ' '^ 



fc. 



a J lap 



(B2) 



By Eq. (jA3l) . ReJCf, ^} « —4. In Appendix [C] it is proved that i/^|-Bi+p = 0{(p'^^^). Consequently, near the critical 
point {hh, fcb), the right hand summand in Ea. (|B2[) is negative so that ai + a2 < 1 as required. 



Appendix C: Complements on the amplitude Ei 

The amplitude of the field £'1+ is a root of the cubic polynomial in Eq. (P5|) which can be solved by Cardano's 
method. Put Y = X + ^ (^) Re{C^}. For the new variable Y, Eq. (|26|) assumes the standard form. 



where, 



Y'^+pY + q = 



p = |i (i«p - 2Remr}) , 9 = II^Mca (4Re{(«)n - mn - ^\Cf 



(CI) 
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As the amplitude Ei is uniquely defined by the system ([6]), it is therefore expected that the cubic in Eq. (jCl[) should 
have a unique real solution in order for the theory to be consistent. The latter holds if and only if the discriminant 

is nonnegative. To prove that D3 > 0, note first that j^^p — 2Re{C^^^} > in the vicinity of a critical point (/if,, fcb). 
This follows from the estimates established in Appendix [K\ Indeed, in the first order of <5o(fcb). 



laefcp - 2Re{(C66)'} 



37r2 



Slik) 



>0 



k—kh 



(C2) 



Next, consider the complex number 



^^27 



2Re{Ce} ( 2Re{(Ce)n " ^l«P ) + *(I«P - ZReKCO^)' 



The positivity condition ()C2p ensures that the coefhcient of the complex number i in the expression of p is indeed 
real. After some algebraic manipulations, it can be shown that. 



D: 



_ ^4 ^ ^ 2 

3 = ^ 



ICP-^-^V 



Thus -D3 > as required. The only real solution Y to Eg. dClj) is then 



Y 




It then follows that. 



\Eu 



V'^+ + ^-' '^± 



^KICP-^^4Re{p}±||CP^^-</'V 



ipl 



■RejCa 



(C3) 



provided t^+t- > 0. The latter condition imposes a limit on the validity of the perturbation theory developed in the 
present study, i.e., the reduction of the system ^ to (|10p is justified if r_|_ + t_ > 0. This is to be expected because 
of the lack of analyticity in Xc of the solution to the nonlinear wave equation ([2]) that can only occur at the critical 
points {hb,kh) at which bound states in the radiation continuum exist. As one gets away from these critical points 
in the (/i, fc)-plane, the solution to the nonlinear wave equation becomes analytic in Xc meaning that all the terms 
that were neglected in finding the principal parts of the amplitudes must now also be taken into account to find a 
solution befitting the series of Eq.Q. The shadowed region depicted in Fig. [2ljb) shows the region of the {h, A;)-plane 
in which the condition t+ + t^ > holds for the first symmetric bound. The presented analysis of the efficiency of 
the second harmonic generation is valid for any choice of the geometrical parameters. Ah = h — h^ and R, and the 
physical parameters, Ec and Xc > 0, which satisfy the conditions (j31a|) and t-|_ + r_ > 0. 
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